Assigned Reading: Geocomputation with R Chapter 2.
Sacramento Mountains, New Mexico
A Visualization of a GIS, by Anne Sexton
The Broad Street Pump, A Cholera Outbreak in London, By: John Snow and digitized by National Geographic
https://www.theguardian.com/news/datablog/2013/mar/15/john-snow-cholera-map
Blue Marble 2012, by Suomi NPP
Ellipsoids
Left: Geoid Cross Section. Right: IGCM Geoid
By International Centre for Global Earth Models (ICGEM) - http://icgem.gfz-potsdam.de/vis3d/longtime / Ince, E. S., Barthelmes, F., Reißland, S., Elger, K., Förste, C., Flechtner, F., Schuh, H. (2019): ICGEM – 15 years of successful collection and distribution of global gravitational models, associated services and future plans. - Earth System Science Data, 11, pp. 647-674,DOI: http://doi.org/10.5194/essd-11-647-2019., CC BY 4.0, https://commons.wikimedia.org/w/index.php?curid=81462823
Reference frame established to represent locations within the frame.
Historically these were locally focused and based on Geoids.
A datum can serve either horizontal (X & Y) or vertical (Z) features.
Components:
Other locations are then measured in relation to the control points.
Number one City Datum. by: Cosmo1976
Control Points Number one City Datum
“The geodetic (or geographic) latitude is the angle between the equatorial plane and the normal (vertical) to the ellipsoid surface at the considered point.” - Proj
Angles on Ellipsoid and Geodetic Coordinates. by: Peter Mercator
By Peter Mercator - Own work, CC BY-SA 3.0, https://commons.wikimedia.org/w/index.php?curid=17717979
\[ \text{Decimal Degrees} = \text{Degrees } + \frac{\text{Minutes}}{60} + \frac{\text{Seconds}}{3600} \] \[ \text{Latitude of Tech} = 42°03'27.7" N \]
\[ 42 + (\frac{03'}{60}) + (\frac{27.7"}{3600}) \]
\[ 42 + 0.05 + 0.007694 = 42.05759 \text{ Decimal Degrees} \]
UTM - Common for field work, planar measurements of meters.
UTM Zones of the Conterminous USA. data by ESRI
Geographic Coordinate Reference System. by: Anna Krystalli
3D to 2d does not work well. But we have worked around this. sorta
Orange Peel. by: Nathan Belz
Many different ways to create two dimensional maps – thousands of projections. None are perfect. There are four main types of projections, each with their own strengths and examples. One of the most common examples of each of these is included in the table below.
Map Projections. by: Daniel Strebe
Unless you really focus on GIS, hardly anything I have said this lecture will matter to you.
Vector data tends to only include features of interest, e.g. bodies of water; whereas a Raster will include, explicitly, the absence of features (e.g. both water and terrestrial areas).
Vector and Raster Data Models
| TAXON | DBH | HEIGHT |
|---|---|---|
| Robinia pseudoacacia | 40 | 24 |
| Quercus alba | 32 | 21 |
all geometries are composed of points.
points only require two coordinates, X & Y.
Y = Latitude (necessary), X = Longitude (necessary)
Z = Elevation (somewhat uncommon)
M = Time or Uncertainty of Measurement (rather uncommon)
SFG is the spatial topology associated with a feature
POINT - A true dimensionless 1 dimensional point
LINESTRING - Sequence of points connected by strings.
POLYGON - Sequence of points connected by strings, which close upon themselves.
i.e. the first and last point of the polygon are the same.
A list column of S3 type containing certain parameters regarding the Geometry/Geometries
Coordinate Reference System (‘crs’)
Precision (‘precision’)
Bounding box (‘bbox’)
Number Empty (‘n_empty’)
Wherein precision refers to the precision of the coordinates forming the geometries. the bbox is a rectangle which includes all pairs of coordinates in the geometry. N_empty, denotes whether the geometry is missing coordinates in a position.
| TAXON | DBH | HEIGHT | LONG | LAT | geometry |
|---|---|---|---|---|---|
| Robinia pseudoacacia | 40 | 24 | -87.67607 | 42.05663 | POINT (-87.67607 42.05663) |
| Quercus alba | 32 | 21 | -87.67399 | 42.05715 | POINT (-87.67399 42.05715) |
Geometry set for 1 feature
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -87.67607 ymin: 42.05663 xmax: -87.67607 ymax: 42.05663
Geodetic CRS: WGS 84
Precision: 2
sfc_POINT of length 1; first list element: 'XY' num [1:2] -87.7 42.1
Released in 2005
First package to hold all major vector types of geometries
Unified Spatial class allowed for support in many spatial statistics packages
Improved mapping of spatial objects
Formed the centerpiece of spatial statistics in R for over a dozen years.
SF holds different geometries in list columns, SP has many different types of objects/classes for holding geometries.
Spatial(Multi)Points
Spatial(Multi)Lines
Spatial(Multi)Polygons
Spatial(Multi)Grids
Less information than SF (‘n_missing’ and ‘precision’ both lacking)
Bounding box
Coordinate Reference System
…Coordinates!
## Formal class 'SpatialPoints' [package "sp"] with 3 slots
## ..@ coords : num [1:15, 1:2] 0.49 0.38 0.27 0.11 0.64 0.03 0.79 0.38 0.28 0.16 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:2] "xc" "yc"
## ..@ bbox : num [1:2, 1:2] 0.03 0.09 0.93 1.78
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "xc" "yc"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr NA
Attachment of attributes, in a data frame, to a Spatial topology can create:
SpatialPointsDataFrame
SpatialLinesDataFrame
SpatialPolgonsDataFrame
SpatialGridsDataFrame
Each S*DF is accessed the same way:
## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 15 obs. of 2 variables:
## .. ..$ temperature: num [1:15] 7.02 4.49 5.4 4.71 4.33 4.73 4.34 5.4 3.92 5.49 ...
## .. ..$ aspect : int [1:15] 16 11 5 17 26 22 38 42 8 35 ...
## ..@ coords.nrs : num(0)
## ..@ coords : num [1:15, 1:2] 0.49 0.38 0.27 0.11 0.64 0.03 0.79 0.38 0.28 0.16 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:2] "xc" "yc"
## ..@ bbox : num [1:2, 1:2] 0.03 0.09 0.93 1.78
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "xc" "yc"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr NA
| temperature | aspect |
|---|---|
| 7.02 | 16 |
| 4.49 | 11 |
| 5.40 | 5 |
| 4.71 | 17 |
| 4.33 | 26 |
| 4.73 | 22 |
We see here that SF contains all of the spatial information in a Simple Feature Collection which can be tidily tucked away in the geometry list column, in SP the analogue to the SFC is the object itself.
SP is mostly geometry with a data slot, SF is mostly data frame with a geometry column.
Western Plant Predictors. Note each cell is the extent of a single raster tile, each tile contains cells.
Note that both cell sizes, and the resolution of the values in cells can, within reason, be converted to finer and coarser resolution. We will discuss these types of calculations next class.
Rasters tend to confuse people. We are going to create our own example here before we talk about them much more.
Fairy shrimp live in bodies of water which dry out in late Spring, and refill in early Fall. We seek to determine where suitable habitat for fairy shrimp are in the Great Basin. In order to do so, we will use satellite imagery to detect valleys which fill with water in Fall, and dry out in late Spring.
Suprise Valley by: John Glen
empty_raster <- raster(
# rasters have 4 bounding edges
# Here we define each 'corner' of the raster'
xmn = 697129.7,
xmx = 811775.7,
ymn = 4388466,
ymx = 4502382,
# Here we set the number of cells 118*118
nrows = 118,
ncols = 118,
# we do NOT manually specify the cell size here.
crs = "+proj=utm +zone=10 +datum=WGS84",
# set the rasters Coordinate Reference System
)
In the code above, we are going to define all four of the essential components of a raster. Clearly, we are defining three explicitly (Bounding Coordinates, Dimensions, CRS), and the remaining one implicitly (Cell resolution).
We can imagine that the raster we are currently creating looks like this. A frame without content. We can see where the bounding coordinates are,
raster_matrix <- matrix(rast_vals_num, # # fill matrix with values,
nrow = empty_raster@nrows, # create matrix same dimensions.
ncol = empty_raster@ncols, # create matrix same dimensions
#ensure filling of matrix goes from upper left to lower right.
byrow = T)
raster_matrix[90:118,112:118] <- 1 # fix the clipping image at edge.
| 0 | 1 | 1 | 0 | 1 | 0 | 0 | 1 |
| 1 | 0 | 1 | 0 | 0 | 1 | 1 | 0 |
| 1 | 0 | 0 | 1 | 1 | 0 | 1 | 0 |
| 0 | 1 | 1 | 0 | 1 | 0 | 0 | 1 |
| 1 | 0 | 1 | 0 | 0 | 1 | 1 | 0 |
| 1 | 0 | 0 | 1 | 1 | 0 | 1 | 0 |
I think it is easiest to imagine that the values within a raster are in the shape of a matrix. Hence, using the code above we could take a vector of values, and bend them, so that each position within the vector matches up with the beginning of a new row.
example_raster_dec <- setValues(empty_raster, raster_matrix)
The width of each raster cell is: 971.57627 meters
The height of each raster cell is: 965.38983 meters
This example_raster_dec contains: 13924 elements
A single raster layer
In general these are rasters which have been classified and we want to extract values from or run calculations with. Now what is great about layers, is that we can do one big thing bricks cannot do, we can load in many layers from many files and create a stack of attributes we are interested in studying on the fly.
The main use of Raster Stacks is to hold layers of similar themes.
Two usage examples:
e.g. yearly mean temperature, mean precipitation, etc.
Do not need be held in memory
As I have mentioned rasters are often generated from satellite imagery; we will discuss rasters which are not developed this way next lecture. Historically most pictures are imaged via the use of three bands. These having spectral values of Red, Green, and Blue (RGB) associated with them. For example our .tif file, is composed of three bands. Note that a Rasterbrick is most often used for loading in these types of imaging data, which can then be processed to form a more typical ‘raster’ dataset.
- do not need to be held in memory,
Focus on sf here
You have scripts to do many types of mapping in your course resources
sf objects are ggplot2 compliant
the order of mapping operations is more important than code
Import example data sets which come with two packages
data("us_states", package = "spData")
us_states <- st_as_sf(us_states)
north_carolina <- read_sf(system.file("shape/nc.shp", package = "sf"))
ggplot(north_carolina) +
geom_sf() +
theme_bw() # we will use this theme for class.
ggplot(north_carolina) +
geom_sf(aes(fill = BIR79)) + # fill the interior of polygons by a variable
theme_bw()
ggplot(north_carolina) +
geom_sf(aes(color = BIR79), # color the borders of polygons by a variable
lwd = 1 # just making the borders thicker.
) +
scale_color_viridis_c(option = "plasma", trans = "sqrt") +
# just using a very colourful color scheme so you can see it.
theme_bw()
ggplot(north_carolina) +
geom_sf(aes(color = BIR79),
fill = NA # set to 'NA' to 'remove' the interior of polygons
) +
scale_color_viridis_c(option = "plasma", trans = "sqrt") +
theme_bw()
ggplot(north_carolina) +
geom_sf(aes(fill = BIR79),
color = NA # set to 'NA' to 'remove' the interior of polygons
) +
scale_color_viridis_c(option = "plasma", trans = "sqrt") +
theme_bw()
ggplot(north_carolina) +
geom_sf(aes(fill = BIR79,
color = FIPS)
) + # both fill and color
guides(color = 'none') + # removed the categorical legend (is too big!)
theme_bw()
ggplot() + # leave empty
geom_sf(data = north_carolina, # put data here
aes(fill = BIR79,
color = FIPS)
) + # both fill and color
guides(color = 'none') + # removed the categorical legend (is too big!)
theme_bw()
ggplot(us_states) + # two data sets.
geom_sf() +
geom_sf(data = north_carolina,
fill = 'purple',
color = NA # hash this line out in your own time to see the difference
) +
theme_bw()
ggplot(north_carolina) + # oh no we drew over North Carolina!
geom_sf(fill = 'purple') +
geom_sf(data = us_states) +
theme_bw()
ggplot() + # two data sets.
geom_sf(data = us_states, fill = NA) +
geom_sf(data = north_carolina, fill = 'purple') +
theme_bw()
ggplot() +
geom_sf(data = us_states) +
geom_sf(data = north_carolina) +
coord_sf(crs = 4267, # we can convert each dataset to the same CRS
datum = NA # we can remove the grid lines/graticules from plot
)
bound <- st_bbox(north_carolina) # retrieve the bbox from the sfc list column
ggplot() +
geom_sf(data = us_states) +
geom_sf(data = north_carolina, aes(fill = BIR79)) +
coord_sf( # use the bbox to 'crop' the extent of the map
xlim = c(bound[1], bound[3]),
ylim = c(bound[2], bound[4])
) +
theme_bw()
For next Lab:
Please install these packages (if you have not done so already):
# install.packages("sf", "raster", "terra", "sp", "tmap", "leaflet", "ggmap")
# optional packages for using parallel processing at a step
# (in our example it won't actually speed up the analyses
# they actually may slow them down!)
# install.packages('snow','parallel')
Download the labs .R script from the course website.
If you are interested in how Drone/Lidar data are collected please check out the ‘RMBL Spatial Data Science Webinar Series’:
For next Lecture: Assigned Reading: Chapter 3 of Spatial Data Science
Future Bonus SDS Office Hours Wednesday Night at 5:00 - 6:00.
Notes on this Lab My lecture notes are in the R script I used to generate all of the novel figures for this presentation. Likewise this presentation is an .HTML file and can be launched from your computer (it was rendered directly from R using the script).
Krystalli, A. https://annakrystalli.me/intro-r-gis/gis.html Accessed 01.20.2022
https://geocompr.robinlovelace.net/spatial-class.html Accessed 01.09.2022
Hijman, R. 05.12.2019 ‘The raster Package’
https://rspatial.org/raster/RasterPackage.pdf Accessed 01.09.2022
https://www.neonscience.org/resources/learning-hub/tutorials/dc-multiband-rasters-r Accessed 01.19.2022
Pebesma, E. https://r-spatial.github.io/sf/articles/sf1.html Accessed 01.10.2022
https://proj.org/operations/conversions/geoc.html Accessed 01.20.2022.
https://rspatial.org/raster/spatial/8-rastermanip.html Accessed 01.11.2022
https://en.wikipedia.org/wiki/Open_Source_Geospatial_Foundation Accessed 01.09.2022
NOAA. What is geodesy? National Ocean Service website, https://oceanservice.noaa.gov/facts/geodesy.html, 1/25/17.
Geocomputation with R. Lovelace, R., Nowosad, J., Muenchow, J. 2022-01-25.
c("raster", "sp", "sf", "tidyverse", "terra") %>%
map(citation) %>%
print(style = "text", na.print = '')
[[1]]
Hijmans R (2022). _raster: Geographic Data Analysis and Modeling_. R
package version 3.5-15, <URL:
https://CRAN.R-project.org/package=raster>.
[[2]]
Pebesma EJ, Bivand RS (2005). "Classes and methods for spatial data in
R." _R News_, *5*(2), 9-13. <URL:
https://CRAN.R-project.org/doc/Rnews/>.
Bivand RS, Pebesma E, Gomez-Rubio V (2013). _Applied spatial data
analysis with R, Second edition_. Springer, NY. <URL:
https://asdar-book.org/>.
[[3]]
Pebesma E (2018). "Simple Features for R: Standardized Support for
Spatial Vector Data." _The R Journal_, *10*(1), 439-446. doi:
10.32614/RJ-2018-009 (URL: https://doi.org/10.32614/RJ-2018-009), <URL:
https://doi.org/10.32614/RJ-2018-009>.
[[4]]
Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R,
Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E,
Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi
K, Vaughan D, Wilke C, Woo K, Yutani H (2019). "Welcome to the
tidyverse." _Journal of Open Source Software_, *4*(43), 1686. doi:
10.21105/joss.01686 (URL: https://doi.org/10.21105/joss.01686).
[[5]]
Hijmans R (2022). _terra: Spatial Data Analysis_. R package version
1.5-18, <URL: https://rspatial.org/terra/>.